function [] = simulateOneExtraProvider(runId)

%%% Add paths
addpath('../model_funcs');
addpath('../model_funcs/dim_specific_funcs/dim3');
addpath('../general_funcs');

%%% Get input / output paths
config = getConfig_consumerAdoptions(runId);

if config.modelId ~= 3
	error('This code is only for modelId == 3 (continuous time models).');
end

outputFolder = createFolderIfNotExist(sprintf('%s/simulation_oneExtraProvider', config.runFolder));
outputFolder1 = createFolderIfNotExist(sprintf('%s/summing_origins', outputFolder));
outputFolder2 = createFolderIfNotExist(sprintf('%s/fixing_destinations', outputFolder));

outputMatFile = sprintf('%s/lambdaDiff.mat', outputFolder);

%%% Load data
load(sprintf('%s/runData.mat', config.runFolder), 'data');
geoLabels    = data.geoLabels;

%%%  Load parameter estimates
load(sprintf('%s/results.mat', config.runFolder), 'params_hat');
params = params_hat; clear params_hat;
% Split params_parts + save everything in object paramsParts
paramsParts  = split_combine_parts(data.dims,  1, params);

%%% Get installed bases + Xdynamic_staticInputs
dynPredInputFile = sprintf('%s/dynamicPredInputs.mat', config.runFolder);
load(dynPredInputFile, 'base');
laggedConsumersOrig      = base.laggedConsumersOrig; % T x K1
laggedConsumersDest      = base.laggedConsumersDest; % T x K2
laggedProviders           = base.laggedProviders;      % T x K2
popAtRisk              = base.popAtRisk';         % K1 x T
Xdynamic_staticInputs  = base.Xdynamic_staticInputs;

%%% Read dimensions
T  = data.dims.dims1{4};
K1 = data.dims.dims1{5};
K2 = data.dims.dims1{6};
NumParams = data.dims.NumParams;
NumParts = length(data.Xparts);
clear base;

% Get spellDurations (and spellWeights to take an average over time)
if isfield(data, 'spellDurations')
	spellWeights = data.spellDurations/sum(data.spellDurations); % T x 1
	spell2period = data.spell2period;
else
	warning('data does not have field spellDurations -- taking uniform weights.');
	spellWeights = 1/T*ones(T,1);
	spell2period = [1:T]';
end

% Get mybeta(coef that multiplies laggedProviders in consumerAdoptions sub-model)
mybeta_idx = find(strcmp(data.Xparts{3}.Xnames, 'log_LaggedProviders_density'));
if isempty(mybeta_idx); error('Could not find beta on log_LaggedProviders_density'); end;
mybeta = paramsParts{3}.beta(mybeta_idx);

% Make object installedBases
laggedConsumersOrig = laggedConsumersOrig'; % K1 x T
laggedConsumersDest = laggedConsumersDest'; % K2 x T
laggedProviders      = laggedProviders';      % K2 x T


%%%%%%%%%%%%%%%%%%%% COMPUTE LAMBDA IN THE BASE CASE (averaged over time) %%%%%%%%%%%%%%%%%%%%
mean_lambdaBase   = zeros(K1, K2);
mean_lambdaCnfact = zeros(K1, K2);

% Assign things that won't depend on spell
data_uu.surfaces                 = data.surfaces;
data_uu.dims.dimType             = data.dims.dimType;
data_uu.dims.NumParts            = data.dims.NumParts;
data_uu.dims.NumParams           = data.dims.NumParams;
data_uu.dims.Xpart_2_NumX        = data.dims.Xpart_2_NumX;
data_uu.dims.Xpart_2_NumX_FEs    = data.dims.Xpart_2_NumX_FEs;
data_uu.dims.Xpart_2_Num_FE_vals = data.dims.Xpart_2_Num_FE_vals;
data_uu.dims.NumFEvals2Keep      = data.dims.NumFEvals2Keep;
data_uu.dims.NumObs   = K1*K2;
data_uu.dims.dims1{1} = K1*K2;
data_uu.dims.dims1{2} = K1;
data_uu.dims.dims1{3} = K2;
data_uu.dims.dims1{4} = 1;
data_uu.dims.dims1{5} = K1;
data_uu.dims.dims1{6} = K2;

for pp = 1:NumParts
	data_uu.Xparts{pp}.NumX_FE_vals = data.Xparts{pp}.NumX_FE_vals;
	data_uu.Xparts{pp}.Xnames       = data.Xparts{pp}.Xnames;
	data_uu.Xparts{pp}.X_FE_labels  = data.Xparts{pp}.X_FE_labels;
end
data_uu.Xparts{1} = data.Xparts{1};
data_uu.Xparts{5} = data.Xparts{5};
data_uu.Xparts{6} = data.Xparts{6};

%%% Modify dimensions of Xdynamic_staticInputs...
%% extraArgs
extraArgs = Xdynamic_staticInputs.extraArgs;
extraArgs = extraArgs{1};
NumPeriods = length(data.periodLabels);
extraArgs = reshape(extraArgs, [NumPeriods K1])';
if ~all(extraArgs - extraArgs(:,1) == 0)
	error('Problem here: extraArgs depends on time...');
end
extraArgs = extraArgs(:,1); % K1 x 1
extraArgs = {extraArgs};
Xdynamic_staticInputs.extraArgs = extraArgs;

%% Xparts_static
Xparts_static = Xdynamic_staticInputs.Xparts_static;

% Xparts{2}.X
X = reshape(Xparts_static{2}.X, [NumPeriods K1 size(Xparts_static{2}.X, 2)]); % NumPeriods x K1 x NumX
X = X(spell2period,:,:); % T x K1 x NumX
X = permute(X, [2 3 1]); % K1 x NumX x T
Xparts_static{2}.X = X;
clear X;

% Xparts{2}.X_FEs
X_FEs = reshape(Xparts_static{2}.X_FEs, [NumPeriods K1 size(Xparts_static{2}.X_FEs, 2)]); % NumPeriods x K1 x NumX_FEs
X_FEs = X_FEs(spell2period,:,:); % T x K1 x NumX_FEs
X_FEs = permute(X_FEs, [2 3 1]); % K1 x NumX_FEs x T
Xparts_static{2}.X_FEs = X_FEs;
clear X_FEs;

% Xparts{3}.X
X = reshape(Xparts_static{3}.X, [NumPeriods K2 size(Xparts_static{3}.X, 2)]); % NumPeriods x K2 x NumX
X = X(spell2period,:,:); % T x K2 x NumX
X = permute(X, [2 3 1]); % K2 x NumX x T
Xparts_static{3}.X = X;
clear X;

% Xparts{3}.X_FEs
X_FEs = reshape(Xparts_static{3}.X_FEs, [NumPeriods K2 size(Xparts_static{3}.X_FEs, 2)]); % NumPeriods x K2 x NumX_FEs
X_FEs = X_FEs(spell2period,:,:); % T x K2 x NumX_FEs
X_FEs = permute(X_FEs, [2 3 1]); % K2 x NumX_FEs x T
Xparts_static{3}.X_FEs = X_FEs;
clear X_FEs;

% Store new result
Xdynamic_staticInputs.Xparts_static = Xparts_static;


% Loop over spells (similar to computeElasticities_ctsTime)
tic
for uu = 1:T
	popAtRisk_uu = popAtRisk(:,uu);       % K1 x 1
	installedBases_uu.laggedConsumersOrig = laggedConsumersOrig(:,uu); % K1 x 1
	installedBases_uu.laggedConsumersDest = laggedConsumersDest(:,uu); % K2 x 1
	data_uu.Xparts{4}.X     = data.Xparts{4}.X(uu,:); % 1 x NumX
	data_uu.Xparts{4}.X_FEs = data.Xparts{4}.X_FEs(uu,:); % 1 x NumX_FEs
	
	Xdynamic_staticInputs_uu.extraArgs     = Xdynamic_staticInputs.extraArgs;
	Xdynamic_staticInputs_uu.Xparts_static = Xdynamic_staticInputs.Xparts_static;
	
	Xdynamic_staticInputs_uu.Xparts_static{2}.X     = Xdynamic_staticInputs.Xparts_static{2}.X(:,:,uu);
	Xdynamic_staticInputs_uu.Xparts_static{2}.X_FEs = Xdynamic_staticInputs.Xparts_static{2}.X_FEs(:,:,uu);
	Xdynamic_staticInputs_uu.Xparts_static{3}.X     = Xdynamic_staticInputs.Xparts_static{3}.X(:,:,uu);
	Xdynamic_staticInputs_uu.Xparts_static{3}.X_FEs = Xdynamic_staticInputs.Xparts_static{3}.X_FEs(:,:,uu);

	
	%%% BASE CASE
	laggedProviders_uu = laggedProviders(:,uu); % K2 x 1
	installedBases_uu.laggedProviders = laggedProviders_uu; % K2 x 1
	
	% Change the Xparts that are "dynamic", based on installed bases that were generated
	Xparts_dynamic_uu = make_Xdynamic_consumerAdoptions(config, Xdynamic_staticInputs_uu, installedBases_uu, true);
	data_uu.Xparts{2} = Xparts_dynamic_uu{2};
	data_uu.Xparts{3} = Xparts_dynamic_uu{3};
	clear Xparts_dynamic_uu;	
	
	% Make logLambdaOffset
	data_uu.logLambdaOffset = log(popAtRisk_uu); % K1 x 1  (this should already be true in the base case)
	data_uu.logLambdaOffset = data_uu.logLambdaOffset - 2e10 * reshape(laggedProviders_uu <= 0, [1 1 K2]); % K1 x 1 x K2
	data_uu.logLambdaOffset = reshape(data_uu.logLambdaOffset, [K1*K2 1]); % (K1*K2) x 1
	
	% Compute lambdaBase
	Vbase_uu = compute_Xbeta2(data_uu.dims, data_uu.Xparts, paramsParts); % (K1*K2) x 1
	logLambdaBase_uu  = data_uu.logLambdaOffset + Vbase_uu;  % (K1*K2) x 1
	clear Vbase_uu;
	lambdaBase_uu     = exp(logLambdaBase_uu); % (K1*K2) x 1
	clear logLambdaBase_uu;
	
	% Increment mean_lambdaBase with results of spell uu
	mean_lambdaBase = mean_lambdaBase + spellWeights(uu) * reshape(lambdaBase_uu,[K1 K2]);
	
	
	%%% CNFACT CASE: add one provider in k2
	laggedProviders_uu = laggedProviders(:,uu) + 1; % K2 x 1
	installedBases_uu.laggedProviders = laggedProviders_uu; % K2 x 1
	
	% Change the Xparts that are "dynamic", based on installed bases that were generated
	Xparts_dynamic_uu = make_Xdynamic_consumerAdoptions(config, Xdynamic_staticInputs_uu, installedBases_uu, true);
	data_uu.Xparts{2} = Xparts_dynamic_uu{2};
	data_uu.Xparts{3} = Xparts_dynamic_uu{3};
	clear Xparts_dynamic_uu;	
	
	% Make logLambdaOffset
	data_uu.logLambdaOffset = log(popAtRisk_uu); % K1 x 1
	data_uu.logLambdaOffset = data_uu.logLambdaOffset - 2e10 * reshape(laggedProviders_uu <= 0, [1 1 K2]); % K1 x 1 x K2
	data_uu.logLambdaOffset = reshape(data_uu.logLambdaOffset, [K1*K2 1]); % (K1*K2) x 1
	
	% Compute lambdaCnfact
	Vcnfact_uu = compute_Xbeta2(data_uu.dims, data_uu.Xparts, paramsParts); % (K1*K2) x 1
	logLambdaCnfact_uu  = data_uu.logLambdaOffset + Vcnfact_uu;  % (K1*K2) x 1
	clear Vcnfact_uu;
	lambdaCnfact_uu     = exp(logLambdaCnfact_uu); % (K1*K2) x 1
	clear logLambdaCnfact_uu;
	
	% Increment mean_lambdaCnfact with results of spell uu
	mean_lambdaCnfact = mean_lambdaCnfact + spellWeights(uu) * reshape(lambdaCnfact_uu,[K1 K2]);

	%%% Display time remaining
	displayRemainingTime(uu, T, 10);
end

%%%%% Compute lambdaDiff
mean_lambdaDiff = mean_lambdaCnfact - mean_lambdaBase; % K1 x K2

if any(mean_lambdaDiff(:) < 0); warning('Something strange: lambda should increase, never decrease!'); end;

%%%%% Save lambdaDiff
save(outputMatFile, 'mean_lambdaBase', 'mean_lambdaCnfact', 'mean_lambdaDiff', '-v7.3');


%%% Compute average total effect of adding a car in k2 across time periods
avgTtlEffect = sum(mean_lambdaDiff,1)'; % K2 x 1

% Display in a table --> normalize by sum to get percentages that is easier to read... (although it's a bit insignificant)
table(geoLabels, 100*avgTtlEffect/sum(avgTtlEffect))

%%% Identify the areas that are the most valuable
[~,idxes] = mink(-avgTtlEffect, 10);

disp('Areas where adding a provider is the most valuable in terms of bringing additional consumers:');
geoLabels(idxes)
%%% Corsica comes super high!!! YES!!! Then Paris


%%%%%% Output that into a table + CSV file
myTable = table(geoLabels, 1e3*avgTtlEffect);
myTable.Properties.VariableNames = {'CantonId', 'avgTtlEffect'};
outputFile1 = sprintf('%s/avgTtlEffect.csv', outputFolder1);
writetable(myTable, outputFile1);
clear myTable;


%%%%%%%%%%% Compute effect of adding a car in Paris (or Lyon / etc.) on consumers in origin cantons k1 %%%%%%%%%%%

paris_idx     = find(strcmp(geoLabels, '75ZZ'));
lyon_idx      = find(strcmp(geoLabels, '69ZZ'));
marseille_idx = find(strcmp(geoLabels, '1398'));
idx_2B01      = find(strcmp(geoLabels, '2B01')); % --> that with the highest value!

my_k2_idxes = [paris_idx, lyon_idx, marseille_idx, idx_2B01];
my_k2_Names = {'Paris', 'Lyon', 'Marseille', 'Canton2B01'};

displayOutFile = sprintf('%s/decomposition.txt', outputFolder2);
fid = fopen(displayOutFile, 'wt');
fclose(fid);

for aa = 1:length(my_k2_idxes)
	myk2 = my_k2_idxes(aa);
	myk2_name = my_k2_Names{aa};
	disp2file(sprintf('%s', myk2_name), displayOutFile, true);

	% Compute effect on origin cantons
	myEffect_avg = mean_lambdaDiff(:,myk2); % K1 x 1
	
	% Compute percentage of the total effect by canton of origin
	pctEffect = 100*myEffect_avg./sum(myEffect_avg);
	disp2file(sprintf('Total effect: %.1f', 1e3*sum(myEffect_avg)), displayOutFile, true);
	disp2file(sprintf('Percentage of total effect that comes from local pop: %.1f percent', pctEffect(myk2)), displayOutFile, true);
	disp2file(sprintf('Percentage of total effect that comes from Paris: %.1f percent', pctEffect(paris_idx)), displayOutFile, true);
	disp2file(sprintf('Percentage of total effect that comes from Bastia: %.4f percent\n', pctEffect(idx_2B01)), displayOutFile, true);

	myTable = table(geoLabels, pctEffect);
	myTable.Properties.VariableNames = {'CantonId', 'pctEffect'};
	outputFile2 = sprintf('%s/%s.csv', outputFolder2, myk2_name);
	writetable(myTable, outputFile2);
	clear myTable;
end

end
